/*
**  File: AquaStrucset.c
**
**  This file is part of the AQUA program suite, which is being
**  developed in the NMR Spectroscopy Department, Bijvoet Center for 
**  Biomolecular Research, Utrecht. 
**
**
**  AUTHOR(S) : Ton Rullmann
**  VERSION   : 1.5.0
**  DATE      : Oct 29, 1997
**
**  This file contains routines for reading, writing and handling
**  (components of) strucsets.
**
**  For safety reasons routines acting on existing strucsets exit
**  if an error occurs, to avoid having a strucset in an ill-defined state.
*/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "Compiler.h"
#include "General.h"

#include "AquaMacros.h"
#include "AquaTypes.h"

/****************  function prototypes  ****************/

#include "AquaFuncts.h"
#include "MenuFuncts.h"
#include "Functs.h"

#include "AquaStrucset.h"


/****************  globals  ****************/

extern Boolean  iaf;                         /* import from Functs */
extern Card     msg_string;                  /* import from Functs */
extern int      terseness;                     /* import from Functs */

static Setlist  *list_head = NULL;
static Card     line, descr, inp;

/* BEGIN_EXTH */

/****************  StrucsetInit  ****************/
/*
**  Initializes a strucset, and associates it with the string <set_name>.
**
**  The address of the strucset is stored in a linked list, headed by <list_head>.
**  The application program should first call StrucsetInit (creating
**  an empty set - emptiness is tested by the macro EMPTYSET; and storing the
**  address for future "call by name"), and then pass the resulting strucset 
**  to other routines, e.g. for reading, selection etc.
*/

void StrucsetInit( Strucset *setptr, char *set_name, const Boolean putlist )
{
    char     *fct_name = "StrucsetInit";
    Strucset *ptr;
    Setlist  *list_new;

    if ( putlist )
    {
	if ( ( ptr = StrucsetFind( set_name ) ) )
	{
	    if ( !EMPTYSET( *ptr ) )
	    {
		printf( " * Will not throw away not-empty set %s\n", set_name );
		exit( 1 );
	    }
	}
	else
	{
	    list_new = CREATE_NEW( Setlist, 1);
	    if ( !list_new )
		NOMEM( "list_new" );
	    list_new->setptr = setptr;
	    list_new->link = list_head;
	    list_head = list_new;
	}
    }

    setptr->molecule = NULL;
    setptr->fst_struc = NULL;
    setptr->lst_struc = NULL;
    setptr->struc_count = 0;
    setptr->atom_count = 0;
    setptr->atom_count_mem = 0;
    setptr->mapped = false;
    setptr->onemap = true;
    setptr->type = undef_data;
    strcpy( setptr->ident, "" );
    strcpy( setptr->derived, "" );
    strcpy( setptr->filnam, "" );
    memcpy( setptr->name, set_name, size_TYPNAM );
    setptr->name[size_TYPNAM-1] = '\0';

    sprintf( msg_string, "Initialization of strucset <%s> completed", setptr->name );
    if ( terseness < 6 ) {
        STATUSm;
    }
}

/****************  WhichStrucset  ****************/
/*
**  Read the name of a strucset and return a pointer to it.
*/

Strucset *WhichStrucset( void )
{
    Card     set_name;
    Menu     *oldmenu;
    Strucset *setptr = NULL;

    if ( list_head->link ) 

/* more than one set defined */
    {
	oldmenu = SetOneMenu( "Give name of strucset to be used (or - to skip):" );
	StrucsetList();
	do
	{
	    MReadWord( set_name );
	    if ( *set_name == '-' || strlen( set_name ) == 0 )
		break;
	    if ( ( setptr = StrucsetFind( set_name ) ) )
		break;
	    printf( " * Strucset <%s> could not be found\n", set_name );
	    MError();
	}
	while 
	    ( true );
	ResetMenu( oldmenu );
    }
    
    else 

/* only one set defined */

	setptr = list_head->setptr;
    
    return( setptr );
}

/****************  StrucsetFind  ****************/
/*
**  Return a pointer to the strucset having name <set_name>.
*/

Strucset *StrucsetFind( char *set_name )
{
    Strucset *setptr;
    Setlist  *list = list_head;
    
    while ( list )
    {
        setptr = list->setptr;
        if ( EQ( setptr->name, set_name ) )
            return( setptr );
        list = list->link;
    }
    return( NULL );
}

/****************  StrucsetList  ****************/
/*
**  List the existing strucsets. All sets are expected to be complete.
*/

void StrucsetList( void )
{
    Strucset *setptr;
    Setlist  *list = list_head;
    size_t   len, max = 0;
    int      atom_count, struc_count;

    if ( terseness < 6 ) {
        puts( "List of defined strucsets:" );
    }
    while ( list )
    {
        len = strlen( list->setptr->name );
        if ( len > max )
            max = len;
        list = list->link;
    }
    list = list_head;
    while ( list )
    {
        setptr = list->setptr;
        atom_count = setptr->atom_count;
        struc_count = setptr->struc_count;
        if ( terseness < 6 ) {
            printf( "<%*s> %s\n", max, setptr->name, setptr->ident );
            printf( "%*s %4i atoms %3i structure(s) %7u bytes\n", max+2, "", 
                    atom_count, struc_count,
                    sizeof(Strucset) + sizeof(Mol) + sizeof(Atom)*(size_t)atom_count + 
                      ( sizeof(Struc) + sizeof(Coord)*(size_t)atom_count )*(size_t)struc_count
                   );
        }
        StrucsetOK( setptr, complete_set );

        list = list->link;
    }
}

/****************  StrucsetStats  ****************/
/*
**  Print some "statistics" of a strucset.
*/

void StrucsetStats( Strucset *setptr )
{
    if ( terseness < 9 ) {
        printf( " * Number of structures: %i\n", setptr->struc_count );
        printf( " * Number of residues:   %i\n", setptr->molecule->res_count );
        printf( " * Number of atoms:      %i\n", setptr->atom_count );
    }
}

/****************  StrucsetCreate  ****************/
/*
**  Allocates memory for a molecule (and possibly structures) in an empty set. 
**  StrucsetInit must have been executed!
**  The set is used properly by:
**   - always referring to set.atom_count
**   - always referring to set.struc_count, or using this loop mechanism:
**     strucptr = set.fst_struc;
**     while ( strucptr ) 
**     {
**         ...
**         strucptr = strucptr->next_struc;
**     }
**  See ReadStrucset, AddAtomToSet, AddCoordToSet for examples.
*/

void StrucsetCreate( Strucset *set, const int atom_count, const int struc_count )
{
    char     *fct_name = "StrucsetCreate";
/*    Mol      *molptr; */
    Struc    *strucptr;
    int      i;

    if ( !EMPTYSET( *set ) )
    {
        printf( " * Will not throw away not-empty set %s\n", set->name );
        exit( 1 );
    }

    set->struc_count = struc_count;
    set->atom_count = atom_count;
    set->molecule = MolCreate( atom_count );
    
    i = 0;
    while ( i++ < struc_count ) 
    {
        strucptr = CREATE_NEW( Struc, 1 );
	if ( !strucptr )
	    NOMEM( "strucptr" );
        if ( i == 1 )
            set->fst_struc = strucptr;
        else
            set->lst_struc->next_struc = strucptr;
        set->lst_struc = strucptr;
        
        strucptr->atom_count = 0;
        strucptr->coords = CREATE_NEW( Coord, atom_count );
        if ( !strucptr->coords )
	    NOMEM( "coords" );
        strucptr->struc_num = i;
    }
}

/****************  MolCreate  ****************/
/*
**  Allocate memory for a molecule (including the atoms).
*/

Mol *MolCreate( const int atom_count )
{
    char  *fct_name = "MolCreate";
    Mol   *molptr;

    molptr = CREATE_NEW( Mol, 1 );    
    if ( !molptr )
	NOMEM( "molecule" );

    molptr->atom_count = 0;                  /* still empty! */
    molptr->atoms = CREATE_NEW( Atom, atom_count );
    if ( !molptr->atoms )
	NOMEM( "atoms" );
    InitAtomS( molptr->atoms, atom_count );

    return( molptr );
}

/****************  MolCopy  ****************/
/*
**  Copy a molecule (not including the atoms and residues).
* 
* Function was a void pointer???
*/

void MolCopy( Mol *new_variable, Mol *old )
{
    Atom     *atomptr;
    Residue  *resptr;

    atomptr = new_variable->atoms;                 /* save ptr */
    resptr  = new_variable->residues;
    memcpy( new_variable, old, sizeof( Mol ) );    /* copy */
    new_variable->atoms = atomptr;                 /* restore ptr */
    new_variable->residues = resptr;
    new_variable->atom_count = 0;                  /* reset atom count */
    
}

/****************  StrucsetCopy  ****************/
/*
**  Copies the contents of strucset <old> to strucset <new>, which should
**  be empty. Memory space in the new set is allocated for <atom_count> atoms
**  by calling StrucsetCreate. <atom_count> can be larger, but not smaller
**  than the number of atoms in the old set.
**  
*/

void StrucsetCopy( Strucset *new_variable, Strucset *old, const int atom_count )
{
    char    *fct_name = "StrucsetCopy";
    Struc   *strucptr, *oldstrucptr, *nextptr;
/*    Atom    *atomptr; */
    Coord   *coordptr;
    
    if ( !StrucsetOK( old, complete_set ) )
        exit( 1 );
    if ( !EMPTYSET( *new_variable ) )
    {
        printf( " * Will not copy into non-empty strucset <%s>\n", new_variable->name );
        exit( 1 );
    }
    if ( atom_count < old->atom_count )
    {
        printf( " * Inconsistent atom_count in copy of strucset <%s> to <%s>\n",
              old->name, new_variable->name );
        printf( "  <%d> in stead of <%d>\n", atom_count, old->atom_count );
        exit( 1 );
    }

    StrucsetCreate( new_variable, atom_count, old->struc_count );
    if ( !StrucsetOK( new_variable, created_set ) )
        exit( 1 );
    
/* copy the Mol (saving the atoms pointer) */

    MolCopy( new_variable->molecule, old->molecule );

/* copy the Struc's (and save the coords and next_struc pointers ) */

    strucptr = new_variable->fst_struc;
    oldstrucptr = old->fst_struc;
    while ( strucptr ) 
    {
        coordptr = strucptr->coords;                         /* save ptrs */
        nextptr = strucptr->next_struc;
        memcpy( strucptr, oldstrucptr, sizeof( Struc ) );    /* copy */
        strucptr->coords = coordptr;                         /* restore ptrs */
        strucptr->next_struc = nextptr;
        strucptr->atom_count = 0;                            /* reset atom count */
        strucptr = nextptr;                                  /* goto next structure */
        oldstrucptr = oldstrucptr->next_struc;
    }

    if ( !StrucsetOK( new_variable, created_set ) )
        exit( 1 );

/* copy the Atom's */

    AddAtomSToSet( new_variable, old->molecule->atoms, old->atom_count );

/* copy the Coord's */

    strucptr = new_variable->fst_struc;
    oldstrucptr = old->fst_struc;
    while ( strucptr ) 
    {
        AddCoordSToSet( new_variable, strucptr, oldstrucptr->coords, old->atom_count );
        strucptr = strucptr->next_struc;
        oldstrucptr = oldstrucptr->next_struc;
    }

    if ( !StrucsetOK( new_variable, part_set ) )
        exit( 1 );

    STATUS( "NEED TO DEFINE RESIDUES AND COPY STRUC_MOL" );
    exit( 1 );
}

/****************  StrucsetJoin  ****************/
/*
**  Joins the contents of strucsets <aa> and <bb> into a new strucset
**  <new>, which should be empty. 
**  Memory space in the new set is allocated by calling StrucsetCreate
**  for a number of atoms that is the sum of the counts of the two existing 
**  sets. The two existing sets should have the same number of structures.
**  At the end MakeResidues is called.
*/

void StrucsetJoin( Strucset *new_variable, Strucset *aa, Strucset *bb )
{
    char    *fct_name = "StrucsetJoin";
    Struc   *strucptr, *oldstrucptr;
    Mol     *molptr;
    
    if ( !StrucsetOK( aa, complete_set ) )
        exit( 1 );
    if ( !StrucsetOK( bb, complete_set ) )
        exit( 1 );
    if ( aa->struc_count != bb->struc_count )
    {
        printf( " * Inconsistent struc_count in joining of strucsets <%s> and <%s>\n",
                aa->name, bb->name );
        printf( "  <%d> in stead of <%d>\n", aa->struc_count, bb->struc_count );
        exit( 1 );
    }

    if ( !EMPTYSET( *new_variable ) )
    {
        printf( " * Will not copy into non-empty strucset <%s>\n", new_variable->name );
        exit( 1 );
    }

    StrucsetCreate( new_variable, aa->atom_count + bb->atom_count, aa->struc_count );
    if ( !StrucsetOK( new_variable, created_set ) )
        exit( 1 );

    if ( !aa->onemap || !bb->onemap )
	ERROR( "IGNORING MODEL DEPENDENT ATOM NAME MAPS" );

/* fill the Mol */

    molptr = new_variable->molecule;
    molptr->inp_filtyp = no_fil_type;

/* fill the Struc's */

/* copy the Atom's */

    AddAtomSToSet( new_variable, aa->molecule->atoms, aa->atom_count );
    AddAtomSToSet( new_variable, bb->molecule->atoms, bb->atom_count );

/* copy the Coord's */

    strucptr = new_variable->fst_struc;
    oldstrucptr = aa->fst_struc;
    while ( strucptr ) 
    {
        AddCoordSToSet( new_variable, strucptr, oldstrucptr->coords, aa->atom_count );
        strucptr = strucptr->next_struc;
        oldstrucptr = oldstrucptr->next_struc;
    }
    strucptr = new_variable->fst_struc;
    oldstrucptr = bb->fst_struc;
    while ( strucptr ) 
    {
        AddCoordSToSet( new_variable, strucptr, oldstrucptr->coords, bb->atom_count );
        strucptr = strucptr->next_struc;
        oldstrucptr = oldstrucptr->next_struc;
    }

    if ( !StrucsetOK( new_variable, complete_set ) )
        exit( 1 );

    MakeResidues( new_variable, true );
}

/****************  StrucsetSort  ****************/
/*
**  Copies the contents of strucset <old> to strucset <new>, which should
**  be empty, and sorts the data of each struc_mol plus corresponding
**  coordinates on the internal atom names as in <sort_mol>. 
**  Effectively it transforms a set for which onemap==false into a new set
**  for which onemap==true.
**  Note that the internal names of the atoms in <sort_mol> must have been
**  set!
**  Memory space in the new set is allocated for the same number of atoms
**  as in the old set.
**
**  In the code <sort_mol> is redefined to the struc_mol of the first
**  model. If we use the strucset molecule some internal atoms will not
**  be found (the strucset molecule has been mapped with a standard AtomLIB, 
**  the first struc_mol has been mapped with a WhatIf-derived AtomLIB; the
**  internal names differ for atoms not occurring in the standard AtomLIB, 
**  such as the terminal NH3 protons). Using the first struc_mol to define
**  the sort order means that any discrepancies (i.e. atom not found) in
**  subsequent models indicate serious inconsistencies between the models.
**  Therefore this generates a fatal error.
*/

void StrucsetSort( Strucset *new_variable, Strucset *old, Mol *sort_mol )
{
    char     *fct_name = "StrucsetSort";
    Struc    *strucptr, *oldstrucptr, *nextptr;
    Coord    *coordptr;
    Mol      *oldmolptr;
    Atom     *at, *fnd;
    int      atom_count = old->atom_count, cnt, idx;
    Boolean  *did_copy, differ = false, sortall = true;

    if ( old->onemap )
	return;

    sort_mol = old->fst_struc->struc_mol;   /* REDEFINE SORT_MOL */

    if ( terseness < 6 ) {
        STATUS( "Starting to sort atoms" );
    }

    if ( !StrucsetOK( old, complete_set ) )
        exit( 1 );
    if ( !EMPTYSET( *new_variable ) )
    {
        printf( " * Will not copy into non-empty strucset <%s>\n", new_variable->name );
        exit( 1 );
    }
    StrucsetCreate( new_variable, atom_count, old->struc_count );
    if ( !StrucsetOK( new_variable, created_set ) )
        exit( 1 );

    did_copy = CREATE_NEW( Boolean, atom_count );
    if ( !did_copy )
	NOMEM( "did_copy" );

/* copy the Mol (saving the atoms pointer) and sort the atoms */
/* atoms are taken from the Mol of the first struc */
/* (since the Mol of the set has not been mapped, cf. MapAtomsInSet !) */

    if ( terseness < 6 ) {
        puts( "Sorting the molecule" );
    }
    MolCopy( new_variable->molecule, old->molecule );
    oldmolptr = old->fst_struc->struc_mol;
    cnt = sort_mol->atom_count;
    at = sort_mol->atoms;
    while ( cnt-- )
    {
	if ( ( fnd = FindAtom1( &idx, at, oldmolptr, use_int, false ) ) )
	{
	    if ( did_copy[idx] )
	    {
		printf( " * Found atom twice: " );
		WriteAtom( stdout, fnd, use_int );
		putchar( '\n' );
		FATAL( "Data error" );
	    }
	    AddAtomToSet( new_variable, fnd );
	    did_copy[idx] = true;
	}
	at++;
    }
/* add remaining atoms (if any) at the end */
    for ( idx=0; idx<atom_count; idx++ )
    {
	if ( !did_copy[idx] )
	{
	    AddAtomToSet( new_variable, &oldmolptr->atoms[idx] );
	    printf( " Atom could not be sorted: " );
	    WriteAtom( stdout, &oldmolptr->atoms[idx], use_int );
	    putchar( '\n' );
	    sortall = false;
	}
	did_copy[idx] = false;
    }

/* copy the Struc's (saving the pointers) and sort the atoms */

    strucptr = new_variable->fst_struc;
    oldstrucptr = old->fst_struc;
    while ( strucptr ) 
    {
        if ( terseness < 6 ) {
	    printf( "Sorting coordinate set %i\n", strucptr->struc_num );
        }
        coordptr = strucptr->coords;                         /* save ptrs */
        nextptr = strucptr->next_struc;
        memcpy( strucptr, oldstrucptr, sizeof( Struc ) );    /* copy */
        strucptr->coords = coordptr;                         /* restore ptrs */
        strucptr->next_struc = nextptr;
        strucptr->atom_count = 0;                            /* reset atom count */
	strucptr->struc_mol = NULL;                          /* no struc_mol */

	oldmolptr = oldstrucptr->struc_mol;
	cnt = sort_mol->atom_count;
	at = sort_mol->atoms;
	while ( cnt-- )
	{
	    if ( (fnd = FindAtom1( &idx, at, oldmolptr, use_int, false ) ))
	    {
		if ( did_copy[idx] )
		{
		    printf( " * Found atom twice: " );
		    WriteAtom( stdout, fnd, use_int );
		    putchar( '\n' );
		    FATAL( "Data error" );
		}
		AddCoordToSet( new_variable, strucptr, &oldstrucptr->coords[idx] );
		did_copy[idx] = true;
		if ( NOTEQID( at->atom_ext, oldmolptr->atoms[idx].atom_ext ) )
		{
		    differ = true;
		    printf( ">>%i ", strucptr->struc_num );
		    PRINTID( at->atom_ext );
		    PRINTID( oldmolptr->atoms[idx].atom_ext );
		    putchar( '\n' );
		}
	    }
	    at++;
	}
/* add remaining atoms (if any) at the end */
/* if sort_mol and struc_mol are inconsistent, the ordering (and therefore
   the name mapping) can be incorrect! */
	for ( idx=0; idx<atom_count; idx++ )
	{
	    if ( !did_copy[idx] )
	    {
		AddCoordToSet( new_variable, strucptr, &oldstrucptr->coords[idx] );
		printf( " Atom could not be sorted: " );
		WriteAtom( stdout, &oldmolptr->atoms[idx], use_int );
		putchar( '\n' );
		sortall = false;
	    }
	    did_copy[idx] = false;
	}
        strucptr = strucptr->next_struc;
        oldstrucptr = oldstrucptr->next_struc;
    }

    if ( !StrucsetOK( new_variable, complete_set ) )
        exit( 1 );

    if ( !sortall ) 
	FATAL( "One or more atoms could not be sorted" );
    if ( !differ )
	puts( " * All models have the same name map (except possibly the unsorted atoms" );
    MakeResidues( new_variable, true );
    new_variable->onemap = true;
}
    
/****************  StrucsetOK  ****************/
/*
**  Check the integrity of a strucset.
*/

Boolean StrucsetOK( Strucset *setptr, const int type )
{
    Mol     *molptr;
    Struc   *strucptr;
    int     atom_count, struc_count = 0;
    
    if ( type == empty_set && !EMPTYSET( *setptr) )
    {
        puts( "set not empty" );
        goto quit;
    }
    if ( type != empty_set && EMPTYSET( *setptr) )
    {
        puts( "set is empty" );
        goto quit;
    }

    molptr = setptr->molecule;
    atom_count = molptr->atom_count;
    
    strucptr = setptr->fst_struc;
    while ( strucptr )
    {
        if ( atom_count != strucptr->atom_count )
        {
            printf( "number of atoms in structure %i inconsistent with molecule definition %i\n",
                    strucptr->atom_count, atom_count );
            goto quit;
        }
        strucptr = strucptr->next_struc;
        struc_count++;
    }
    if ( struc_count != setptr->struc_count )
    {
        printf( "number of structures %i inconsistent with set definition %i\n",
              struc_count, setptr->struc_count );
        goto quit; 
    }

    if ( atom_count != setptr->atom_count )
        if (    !( type == created_set && atom_count == 0 ) 
             && !( type == part_set ) )
        {
            printf( "number of atoms in molecule %i inconsistent with set definition %i\n",
                    atom_count, setptr->atom_count );
            goto quit;
        }
    
    strucptr = setptr->fst_struc;
    while ( strucptr )
    {
        if ( strucptr->atom_count != setptr->atom_count )
            if (    !( type == created_set && strucptr->atom_count == 0 ) 
                 && !( type == part_set ) )
            {
                printf( "number of atoms in structure %i inconsistent with set definition %i\n",
                        strucptr->atom_count, setptr->atom_count );
                goto quit;
            }
        strucptr = strucptr->next_struc;
    }

    return( true );

quit:
    printf( " * ERROR(S) for set <%s>\n", setptr->name );
    return( false ); 
}

/****************  DefaultStrucsetIdent  ****************/
/*
**  Set the ident field of a Strucset to the default value
**  determined from the type field.
*/

void DefaultStrucsetIdent( Strucset *set )
{
    memcpy( set->ident, QtypTxt( set->type ), size_LINE );
    set->ident[size_LINE-1] = '\0';
}

/****************  AddAtomToSet  ****************/
/*
**  
*/

void AddAtomToSet( Strucset *setptr, Atom *atomptr )
{
    Mol *molptr;
    
    molptr = setptr->molecule;
    if ( molptr->atom_count >= setptr->atom_count )
    {
        printf( " * Too many atoms in molecule (current: %i, new: %i, set: %i)\n",
                molptr->atom_count, molptr->atom_count + 1, setptr->atom_count );
        exit( 1 );
    }
    memcpy( &molptr->atoms[molptr->atom_count++], atomptr, sizeof(Atom) );
}

/****************  AddAtomSToSet  ****************/
/*
**  
*/

void AddAtomSToSet( Strucset *setptr, Atom *atomptr, const int count )
{
/*    char *fct_name = "AddAtomSToSet"; */
    Mol  *molptr;
    
    molptr = setptr->molecule;
    if ( molptr->atom_count + count > setptr->atom_count )
    {
        printf( " * Too many atoms in molecule (current: %i, new: %i, set: %i)\n",
                molptr->atom_count, molptr->atom_count + count, setptr->atom_count );
        exit( 1 );
    }
/*
    sprintf( msg_string, "copying atoms %i through %i into molecule",
             molptr->atom_count + 1, molptr->atom_count + count );
    STATUSm;
*/
    memcpy( &molptr->atoms[molptr->atom_count], atomptr, (size_t)count*sizeof(Atom) );
    molptr->atom_count += count;
}

/****************  AddCoordToSet  ****************/
/*
**  
*/

void AddCoordToSet( Strucset *setptr, Struc *strucptr, Coord *coordptr )
{
/*    char *fct_name = "AddCoordToSet"; */
    
    if ( strucptr->atom_count >= setptr->atom_count )
    {
        printf( " * Too many atoms in structure (current: %i, new: %i, set: %i)\n",
                strucptr->atom_count, strucptr->atom_count + 1, setptr->atom_count );
        exit( 1 );
    }
/*
    sprintf( msg_string, "copying atom %i into structure", strucptr->atom_count + 1 );
    STATUSm;
*/
    memcpy( &strucptr->coords[strucptr->atom_count++], coordptr, sizeof(Coord) );
}

/****************  AddCoordSToSet  ****************/
/*
**  
*/

void AddCoordSToSet( Strucset *setptr, Struc *strucptr, Coord *coordptr, 
                     const int count )
{
/*    char *fct_name = "AddCoordSToSet"; */
    
    if ( strucptr->atom_count + count > setptr->atom_count )
    {
        printf( " * Too many atoms in structure (current: %i, new: %i, set: %i)\n",
                strucptr->atom_count, strucptr->atom_count + count, setptr->atom_count );
        exit( 1 );
    }
/*
    sprintf( msg_string, "copying atoms %i through %i into structure",
             strucptr->atom_count + 1, strucptr->atom_count + count );
    STATUSm;
*/
    memcpy( &strucptr->coords[strucptr->atom_count], coordptr, (size_t)count*sizeof(Coord) );
    strucptr->atom_count += count;
}

/****************  AddAtomSToMol  ****************/
/*
**  
*/

void AddAtomSToMol( Mol *molptr, Atom *atomptr, const int count, const int set_atom_count )
{
/*    char *fct_name = "AddAtomSToMol"; */
    
    if ( molptr->atom_count + count > set_atom_count )
    {
        printf( " * Too many atoms in molecule (current: %i, new: %i, set: %i)\n",
                molptr->atom_count, molptr->atom_count + count, set_atom_count );
        exit( 1 );
    }
/*
    sprintf( msg_string, "copying atoms %i through %i into molecule",
             molptr->atom_count + 1, molptr->atom_count + count );
    STATUSm;
*/
    memcpy( &molptr->atoms[molptr->atom_count], atomptr, (size_t)count*sizeof(Atom) );
    molptr->atom_count += count;
}

/****************  DeleteMol  ****************/
/*
*/

void DeleteMol( Mol *molptr )
{
    if ( molptr->residues )
    {
	free( molptr->residues->atomidxs );
	free( molptr->residues->atomptrs );
	free( molptr->residues );
    }
    free( molptr->atoms );
    free( molptr );
}

/****************  ProcessStrucsetFile  ****************/
/*
**  Read a strucset file for which the name is known
**  (simplified version of ReadAquaSet)
*/

Boolean ProcessStrucsetFile( Strucset *set, const int strucs_to_alloc )
{
    FILE     *fp;
    Boolean  ok = false;
     
    if ( terseness < 6 ) {
        printf( "Now opening: %s\n", set->filnam );
    }
    if ( NOTOPENr( fp, set->filnam ) )
    {
        printf( " * Strucset-file <%s> could not be opened\n", set->filnam );
        if ( !iaf )
	    exit( 1 );
	goto end;
    }

    if ( !ReadAquaIDLines( fp ) )
    {
	fclose( fp );
	goto end;
    }
    ok = ReadStrucset( fp, set, strucs_to_alloc );
    fclose( fp );

end:   
    return( ok );
}

/****************  ReadAquaSet  ****************/
/*
**  Open an Aqua set file, and call ReadStrucset.
*/

Boolean ReadAquaSet( Strucset *set, const int strucs_to_alloc )
{
    FILE     *fp;
    Qtyp     qtype;
    Filnam   filnam;
    Menu     *oldmenu;
    Boolean  ok = false;
    int      item;
    
    static MenuItem menulist[] =
    {
        { dummy_str_data, "dummy", none, on },
        { orig_str_data,  "orig",  none, on }, 
        { explo_str_data, "explo", none, on }, 
        { pseu_str_data,  "pseu",  none, on }, 
        { sel_str_data,   "sel",   none, off }, 
        { item_END,       "end",  none, on } 
    };
    static Menu menu =
    {
        "ReadAquaSet menu", "file selector", "Give selector:", 
        sizeof(menulist)/sizeof(MenuItem)
    };
    
    oldmenu = SetMenu( &menu, menulist );

    while ( true )
    {
        item = ReadMenuItem();
        switch ( item )
        {
            case dummy_str_data:
            case orig_str_data:
            case explo_str_data:
            case pseu_str_data:
            case sel_str_data:
                qtype = (Qtyp) item;
                goto data;
            case item_END:
            case item_NULL:
                goto end;
            default:
                ProcessMenuError( menu_fatal );
                break;
        }
    }
data:
    fp = OpenAquaFile( qtype, "r", filnam, false );
    if ( !fp )
	goto end;

    if ( !ReadAquaIDLines( fp ) )
    {
	fclose( fp );
	goto end;
    }
    ok = ReadStrucset( fp, set, strucs_to_alloc );
    set->type = qtype;
    strcpy( set->filnam, filnam );
    if ( IsEmptyStr( set->ident ) )
	DefaultStrucsetIdent( set );
    fclose( fp );

end:
    ResetMenu( oldmenu );
    return( ok );
}

/****************  WriteAquaSet  ****************/
/*
**  Open an Aqua set file, and call WriteStrucset.
*/

Boolean WriteAquaSet( Strucset *set, const Boolean toscreen, const Boolean reuse )
{
    FILE    *fp;
    Qtyp    qtype;
    Filnam  filnam;
    
    qtype = set->type;
    if ( toscreen )
    {
	WriteStrucset( stdout, set );
	return( true );
    }

    fp = OpenAquaFile( qtype, "w", filnam, reuse );
    if ( fp )
    {
        WriteAquaIDLines( fp, qtype, filnam );
        WriteStrucset( fp, set );
        strcpy( set->filnam, filnam );
	fclose( fp );
	return( true );
    }
    else
	return( false );
}

/****************  ReadStrucset  ****************/
/*
**  Process the file header, allocate memory (if necessary),
**  and read the molecule (ReadMol) and the structures (ReadStruc).
**
**  Space is allocated for <strucs_to_alloc> structures. 
**  If <strucs_to_alloc> is negative, 
**  the struc_count value on the file is used instead.
**
**  At the end MakeResidues is called.
*/

Boolean ReadStrucset( FILE *fp, Strucset *set, int strucs_to_alloc )
{
    Struc    *strucptr;
    Boolean  done = false, ok = false;
    Card     ident;
    int      count = 0, atom_count = 0, struc_count = 0;
        
/* process file header */

    while ( !done )
    {
	if ( ! ProcessAquaLine( line, descr, inp, size_LINE, fp ) )
	{
	    if ( !iaf )
		exit( 1 );
	    goto end;
	}
        if      ( EQn( descr, "IDENT", 5 ) )
            strcpy( ident, inp );
        else if ( EQn( descr, "ATOM_COUNT", 10 ) )
            sscanf( inp, "%i", &atom_count );
        else if ( EQn( descr, "STRUC_COUNT", 11 ) )
            sscanf( inp, "%i", &struc_count );
        else if ( EQn( descr, "MOLECULE", 8 ) )      /* start of molecule */
            done = true;
    }

    if ( !atom_count )
    {
        puts( " * Data missing from file header" );
        if ( !iaf )
	    exit( 1 );
	goto end;
    }

/* replace a negative strucs_to_alloc value by struc_count */

    if ( strucs_to_alloc < 0 )
	strucs_to_alloc = struc_count;

    if ( strucs_to_alloc > struc_count )
    {
        puts( " * Additional strucset input not yet implemented" );
        if ( !iaf )
	    exit( 1 );
	goto end;
    }


    if ( EMPTYSET( *set ) )
    {
/* create the set, i.e. allocate memory, set all pointers etc. */
/* fill the molecule - data from first input file */

        StrucsetCreate( set, atom_count, strucs_to_alloc );
        strcpy( set->ident, ident );
        ReadMol( fp, set->molecule, atom_count );
	set->mapped = true;
	MakeResidues( set, true );
    }
    else
    {
/* check the available memory space */
/* compare stored molecule with data from next input file */

        puts( " * Additional strucset input not yet implemented" );
        if ( !iaf )
	    exit( 1 );
	goto end;
    }
    if ( terseness < 6 ) {
        puts( "A molecule has been read" );
    }


/* fill the new structures (as many as have been allocated; the file
   may contain more) */

    strucptr = set->fst_struc;
    while ( strucptr ) 
    {
        strucptr = ReadStruc( fp, strucptr, atom_count );
        count++;
        if ( terseness < 6 ) {
	    printf( "Coordinate set %i has been read\n", count );
        }
    }

    if ( terseness < 6 ) {
        printf( "%i Coordinate set%s been read\n", count, 
            count == 1 ? " has" : "s have" );
    }
    ok = true;

end:    
    return( ok );
}

/****************  WriteStrucset  ****************/
/*
**  
*/

void WriteStrucset( FILE *fp, Strucset *set )
{
    Struc   *strucptr;
        
    fprintf( fp, "$ IDENT: %s\n", set->ident );
    fprintf( fp, "$ ATOM_COUNT: %i\n", set->atom_count );
    fprintf( fp, "$ STRUC_COUNT: %i\n", set->struc_count );

    if ( EMPTYSET( *set ) )
        return;

    WriteMol( fp, set->molecule );

    strucptr = set->fst_struc;
    while ( strucptr )
        strucptr = WriteStruc( fp, strucptr );
    if ( terseness < 6 ) {
        printf( "%i structure%s been written\n", set->struc_count, 
                set->struc_count == 1 ? " has" : "s have" );
    }
}

/****************  ReadMol  ****************/
/*
**  
*/

void ReadMol( FILE *fp, Mol *mol, const int count )
{
    Boolean done = false;
    int temp;
    
    while ( !done )
    {
	if ( ! ProcessAquaLine( line, descr, inp, size_LINE, fp ) )
	    exit( 1 );
        if      ( EQn( descr, "From file", 9 ) )
            sscanf( inp, format_FILNAM, mol->inp_filnam );
        else if ( EQn( descr, "File type", 9 ) ) {
            /* Just to be explicit for comiler warnings */
            sscanf( inp, "%i", &temp );
            mol->inp_filtyp = (Filtyp) temp;
        } else if ( EQn( descr, "Ident", 5 ) )
            strcpy( mol->ident, inp );
        else if ( EQn( descr, "Atom count", 10 ) )
        {
            sscanf( inp, "%i", &mol->atom_count );
            if ( mol->atom_count != count )
            {
                puts( " * Inconsistent atom_count in molecule" );
                printf( "  <%d> in stead of <%d>\n",
		 mol->atom_count, count );
                exit( 1 );
            }
        }
        else if ( EQn( descr, "BEGIN_DATA", 10 ) )
            done = true;
    }
    ReadAtoms( fp, mol->atoms, mol->atom_count );
}

/****************  WriteMol  ****************/
/*
**  
*/

void WriteMol( FILE *fp, Mol *mol )
{
    fprintf( fp, "$ MOLECULE:\n" );
    fprintf( fp, "$ From file: %s\n", mol->inp_filnam );
    fprintf( fp, "$ File type: %i (%s)\n", mol->inp_filtyp, FiltypNam( mol->inp_filtyp ) );
    fprintf( fp, "$ Ident: %s\n", mol->ident );
    fprintf( fp, "$ Atom count: %i\n", mol->atom_count );
    fprintf( fp, "$ BEGIN_DATA: \n" );

    WriteAtoms( fp, mol->atoms, mol->atom_count );
}

/****************  ReadAtoms  ****************/
/*
**  
*/

void ReadAtoms( FILE *fp, Atom *atoms, const int atom_count )
{
    Status  stat;
    int     count = 0;

    while ( ++count <= atom_count )
    {
        stat = ReadAtom( fp, atoms++ );
	if ( stat != status_ok )
	{
	    printf( " * Error occurred after reading %i atoms\n", count - 1 );
	    exit( 1 );
	}
    }
}

/****************  WriteFAtoms  ****************/
/*
**  
*/

void WriteFAtoms( FILE *fp, Atom *atoms, int count )
{
    while ( count-- )
    {
        WriteFAtom( fp, atoms++, use_both );
        fputc( '\n', fp );
    }
}

/****************  WriteAtoms  ****************/
/*
**  
*/

void WriteAtoms( FILE *fp, Atom *atoms, int count )
{
    while ( count-- )
    {
        WriteAtom( fp, atoms++, use_both );
        fputc( '\n', fp );
    }
}

/****************  ReadAtom  ****************/
/*
**  
*/

Status ReadAtom( FILE *fp, Atom *atom )
{
    Card     line;
    Status   stat;

    if ( !ReadWord( line, size_LINE, false, fp ) )  /* input first word */
	return( eof_err );

    if ( *line == CHAIN )
    {
	CopyChars( atom->chain, line+1, size_ID );      /* read chain name */

	if ( !ReadWord( line, size_LINE, false, fp ) )  /* input next word */
	{
	    stat = eof_err;
	    goto err;
	}
    }

    CopyChars( atom->res_nam, line, size_ID );      /* read residue name */

    if ( fscanf( fp, "%i%c", &atom->res_num_orig, &atom->res_insert ) < 2 )
                                                    /* read orig residue number */
                                                    /* and insertion code */
    {
        puts( " * Error in reading residue number from Aqua file" );
        stat = read_err;
	goto err;
    }
/*--atom->res_num = atom->res_num_orig;  */         /* res_num may be changed later.... */

    if ( !ReadWord( line, size_LINE, false, fp ) )  /* input next word */
    {
	stat = eof_err;
	goto err;
    }

    CopyChars( atom->atom_ext, line, size_ID );     /* read ext atom name */
    
    if ( !ReadWord( line, size_LINE, false, fp ) )  /* input next word */
    {
	stat = eof_err;
	goto err;
    }

    CopyChars( atom->atom_int, line, size_ID );     /* read int atom name */
    
    return( status_ok );
    
err:
    printf( " * Last read: %s\n", line );
    return( stat );
}

/****************  WriteFAtom  ****************/
/*
**  Write atom identification in a fixed length format.
*/

void WriteFAtom( FILE *fp, Atom *atom, const int which )
{
    size_t  l;

    fputc( CHAIN, fp );
    if ( ( l = LenChars( atom->chain, size_ID ) ) )
    {
	FPrintChars( fp, atom->chain, l );
        fputc( ' ', fp );
    }
    else
        fprintf( fp, "  " );
    if ( fp == stdout )
        fprintf( fp, "(%3i) ", atom->res_num_new );
    FPRINTID( fp, atom->res_nam );
    fprintf( fp, " %3i%c ", atom->res_num_orig, atom->res_insert );
    if ( which != use_int )
	FPRINTID( fp, atom->atom_ext );
    if ( which == use_both )
	fputc( ' ', fp );
    if ( which != use_ext )
	FPRINTID( fp, atom->atom_int );
}

/****************  WriteAtom  ****************/
/*
**  Write atom identification in a variable length format.
*/

void WriteAtom( FILE *fp, Atom *atom, const int which )
{
    size_t  l;

    if ( ( l = LenChars( atom->chain, size_ID ) ) )
    {
        fputc( CHAIN, fp );
	FPrintChars( fp, atom->chain, l );
        fputc( ' ', fp );
    }
    if ( fp == stdout )
        fprintf( fp, "(%i) ", atom->res_num_new );
    FPRINTIDv( fp, atom->res_nam );
    fprintf( fp, " %i%c", atom->res_num_orig, atom->res_insert );
    if ( atom->res_insert != ' ' )
	fputc( ' ', fp );
    if ( which != use_int )
	FPRINTIDv( fp, atom->atom_ext );
    if ( which == use_both )
	fputc( ' ', fp );
    if ( which != use_ext )
	FPRINTIDv( fp, atom->atom_int );
}

/****************  WriteFResidue  ****************/
/*
**  Write residue identification in a fixed length format
*/

void WriteFResidue( FILE *fp, Residue *residue )
{
    size_t  l;

    fputc( CHAIN, fp );
    if ( (l = LenChars( residue->chain, size_ID ) ))
    {
	FPrintChars( fp, residue->chain, l );
        fputc( ' ', fp );
    }
    else
        fprintf( fp, "  " );
    if ( fp == stdout )
        fprintf( fp, "(%3i) ", residue->res_num_new );
    FPRINTID( fp, residue->res_nam );
    fprintf( fp, " %3i%c ", residue->res_num_orig, residue->res_insert );
}

/****************  WriteResidue  ****************/
/*
**  Write residue identification in a variable length format
*/

void WriteResidue( FILE *fp, Residue *residue )
{
    size_t  l;

    if ( ( l = LenChars( residue->chain, size_ID ) ) )
    {
        fputc( CHAIN, fp );
	FPrintChars( fp, residue->chain, l );
        fputc( ' ', fp );
    }
    if ( fp == stdout )
        fprintf( fp, "(%i) ", residue->res_num_new );
    FPRINTIDv( fp, residue->res_nam );
    fprintf( fp, " %i%c", residue->res_num_orig, residue->res_insert );
    if ( residue->res_insert != ' ' )
	fputc( ' ', fp );
}

/****************  ReadStruc  ****************/
/*
**  
*/

Struc *ReadStruc( FILE *fp, Struc *struc, const int count )
{
    Boolean done = false;
    
    while ( !done )
    {
	if ( ! ProcessAquaLine( line, descr, inp, size_LINE, fp ) )
	    exit( 1 );
        if      ( EQn( descr, "STRUCTURE", 9 ) )
            sscanf( inp, "%i", &struc->struc_num );
        else if ( EQn( descr, "From file", 9 ) )
            sscanf( inp, format_FILNAM, struc->inp_filnam );
        else if ( EQn( descr, "File structure number", 21 ) )
            sscanf( inp, "%i", &struc->fil_struc_num );
        else if ( EQn( descr, "Ident", 5 ) )
            strcpy( struc->ident, inp );
        else if ( EQn( descr, "Atom count", 10 ) )
        {
            sscanf( inp, "%i", &struc->atom_count );
            if ( struc->atom_count != count )
            {
                puts( " * Inconsistent atom_count in structure" );
                printf( "  <%d> in stead of <%d>\n",
		 struc->atom_count, count );
                exit( 1 );
            }
        }
        else if ( EQn( descr, "BEGIN_DATA", 10 ) )
            done = true;
    }
    
    ReadCoords( fp, struc->coords, struc->atom_count );

    return( struc->next_struc );
}

/****************  WriteStruc  ****************/
/*
**  
*/

Struc *WriteStruc( FILE *fp, Struc *struc )
{
    fprintf( fp, "$ STRUCTURE: %i\n", struc->struc_num );
    fprintf( fp, "$ From file: %s \n", struc->inp_filnam);
    fprintf( fp, "$ File structure number: %i\n", struc->fil_struc_num );
    fprintf( fp, "$ Ident: %s\n", struc->ident );
    fprintf( fp, "$ Atom count: %i\n", struc->atom_count );
    fprintf( fp, "$ BEGIN_DATA: \n" );

    WriteCoords( fp, struc->coords, struc->atom_count );
    
    return( struc->next_struc );
}

/****************  ReadCoords  ****************/
/*
**  
*/

void ReadCoords( FILE *fp, Coord *coords, const int atom_count )
{
    int  count = 0;

    while ( ++count <= atom_count )
    {
        if ( fscanf( fp, "%8f%8f%8f\n", &coords->x, &coords->y, &coords->z )
	     < 3 )
	{
	    printf( " * Error occurred after reading coordinates of %i atoms\n", 
	            count - 1 );
	    exit( 1 );
	}
        coords++;
    }
}

/****************  WriteCoords  ****************/
/*
**  
*/

void WriteCoords( FILE *fp, Coord *coords, int count )
{
    while ( count-- )
    {
        fprintf( fp, "%8.3f%8.3f%8.3f\n", coords->x, coords->y, coords->z );
        coords++;
    }
}

/****************  RenumResidues  ****************/
/*
**  Renumber residues by copying the new, sequential residue number
**  to the original one, and blanking the chain identifier and 
**  residue insertion code.
*/

void RenumResidues( Strucset *set )
{
    char    *fct_name = "RenumResidues";
    Mol     *molecule = set->molecule;
    Atom    *atoms = molecule->atoms;
    int     count = molecule->atom_count;
    
    if ( ! molecule->residues )
	FATAL( "no residues defined yet" );

    while ( count-- )
    {
	atoms->res_num_orig = atoms->res_num_new;
	atoms->res_insert = ' ';
	Blank( atoms->chain, size_ID );
	atoms++;
    }
    free( molecule->residues );
    MakeResidues( set, true );
}
